Monotone TDLNM
Daniel Mork, Harvard T.H. Chan School of Public Health
Aug 15, 2023
Create lagged data
We create a lagged dataset with current exposures and exposures from 20 days prior. We differentiate the dataset to model summer (May-Sept) and winter (Oct-Apr) separately because very cold and very warm temperatures are expected to both have detrimental effects on health and we wish to model a monotone effect.
As death counts are large, we assume log-death rate to be a continuous variable under a standard linear model framework. We control for time using an interaction between year and month as well as day of week. We compare three methods: Monotone-TDLNM, TDLNM, and GAM under a penalized spline DLNM framework.
data("chicagoNMMAPS")
setDT(chicagoNMMAPS)
lags <- 20
chicagoNMMAPS[, paste0("temp.", 0:lags) := shift(temp, 0:lags, type = "lag")]
cols <- c("death", "month", "year", "dow", paste0("temp.", 0:lags))
data <- chicagoNMMAPS[complete.cases(chicagoNMMAPS[, ..cols])]
data$logdeath <- log(data$death)
par(mfrow = c(1, 2))
hist(data$death, xlab = "Daily deaths", main = "")
hist(data$logdeath, xlab = "Log daily deaths", main = "")
cols <- paste0("temp.", 0:lags)
# Summer data
summer_data <- data[month > 4 & month < 10]
summer_temp_dat <- as.matrix(summer_data[, ..cols])
temp_splits <- seq(min(summer_temp_dat), max(summer_temp_dat), length.out = 20)
# N days: summer
summer_data[, .N]## [1] 2142# Winter data
winter_data <- data[month < 5 | month > 9]
winter_temp_dat <- -as.matrix(winter_data[, ..cols])
w_temp_splits <- seq(min(winter_temp_dat), max(winter_temp_dat), length.out = 20)
# N days: winter
winter_data[, .N]## [1] 2952Summer model
We model summer mortality and lagged temperature data using observations from May through September. For Monotone-TDLNM we use informative priors to incorporate knowledge from past research indicating high temperatures from up to 5 days prior may be related to increased mortality. Specifically, we assign a prior for lags 0-5 that we have 95% confidence the probability of an effect lies between 0.8 and 0.99. For the remaining lags we indicate with 95% confidence that the probability of an effect is between 0.01 and 0.8.
Monotone-TDLNM
# g0 <- rep(0, lags + 1) # prior prob 0.025-0.975
# g0[1:6] <- 2.991 # prior prob 0.8-0.99
# s0 <- diag(lags + 1) * 1.869^2
# for (s in 1:6)
# s0[s, s] <- 0.819^2
g0 <- rep(0, lags + 1) # prior prob 0.005-0.995
g0[1:6] <- 4.119 # prior prob 0.8-0.99
s0 <- diag(lags + 1) * 2.701^2
for (s in 1:6)
s0[s, s] <- 0.599^2
ts0 <- rep(1, lags)
ts0[1:6] <- 10
ts0 <- ts0 / sum(ts0)
mlist <-
foreach(m = 1:restarts, .errorhandling = "remove",
.packages = "dlmtree",
.verbose = FALSE) %dopar%
{
set.seed(m * 100)
tdlnm(logdeath ~ as.factor(year) * as.factor(month) + as.factor(dow),
data = summer_data,
exposure.data = summer_temp_dat,
exposure.splits = temp_splits,
n.trees = trees, n.burn = burn, n.iter = iter, n.thin = 10,
monotone = T, verbose = F, altmcmc = 1,
monotone.gamma0 = g0, monotone.sigma = s0, time.split.prob = ts0
)
}
summer_deaths <- combine.models(mlist)
s_summer_deaths <- summary(summer_deaths, cenval = 20, pred.at = 10:32,
mcmc = TRUE, verbose = F)
s_summer_deaths$matfit <- out_trans(s_summer_deaths$matfit)
s_summer_deaths$cilower <- out_trans(s_summer_deaths$cilower)
s_summer_deaths$ciupper <- out_trans(s_summer_deaths$ciupper)
s_summer_deaths$plot.dat$Est <- out_trans(s_summer_deaths$plot.dat$Est)
# Posterior prob of nonzero effect and time splits
par(mfrow = c(1, 2))
boxplot(1/(1+exp(-summer_deaths$zirtGamma)), ylim = c(0, 1), main = "Prob. nonzero ERC", names = 0:lags, cex = 0.5); points(1/(1+exp(-g0-1.96*sqrt(diag(s0)))), col = "red", pch = "-", cex = 3); points(1/(1+exp(-g0+1.96*sqrt(diag(s0)))), col = "red", pch = "-", cex = 3);
boxplot(summer_deaths$timeProbs, ylim = c(0, 1), main = "Prob. time split", names = sapply(1:lags, function(i) paste0(i-1, "/", i)), cex = 0.5); points(ts0, col = "blue", pch = 16)
# Linear model assumptions
par(mfrow = c(1, 2))
plot(summer_deaths$Yhat,
summer_data$logdeath - summer_deaths$Yhat,
xlab = "Fitted", ylab = "Residuals"); abline(h = 0)
qqnorm(scale(summer_data$logdeath - summer_deaths$Yhat, center = F), main = ""); abline(0, 1)
# Convergence
s <- lapply(mlist, summary, cenval = -20, pred.at = 10:32, mcmc = TRUE, verbose = FALSE)
summer_mcmc_mtdlnm <- do.call(mcmc.list, lapply(1:restarts, function(i) {
d <- as.data.frame.table(s[[i]]$dlm_mcmc)
d$Var1 <- paste0("Lag", rep(rep(0:lags, times = length(s[[i]]$pred.at)), times = iter/10))
d$Var2 <- paste0("Temp", rep(rep(s[[i]]$pred.at, each = length(0:lags)), times = iter/10))
setDT(d)
mcmc(dcast(d[Var2 != "Temp20"], Var3 ~ Var1 + Var2, value.var = "Freq")[, Var3 := NULL][])
}))
gr_summer_mtdlnm <- gelman.diag(summer_mcmc_mtdlnm, transform = F, autoburnin = F, multivariate = F)TDLNM
mlist <-
foreach(m = 1:restarts, .errorhandling = "remove",
.packages = "dlmtree",
.verbose = FALSE) %dopar%
{
set.seed(m * 100)
tdlnm(logdeath ~ as.factor(year) * as.factor(month) + as.factor(dow),
data = summer_data,
exposure.data = summer_temp_dat,
exposure.splits = temp_splits,
n.burn = burn, n.iter = iter, n.thin = 10,
verbose = F)
}
summer_deaths2 <- combine.models(mlist)
s_summer_deaths2 <- summary(summer_deaths2, cenval = 20, pred.at = 10:32,
mcmc = TRUE, verbose = F)
s_summer_deaths2$matfit <- out_trans(s_summer_deaths2$matfit)
s_summer_deaths2$cilower <- out_trans(s_summer_deaths2$cilower)
s_summer_deaths2$ciupper <- out_trans(s_summer_deaths2$ciupper)
s_summer_deaths2$plot.dat$Est <- out_trans(s_summer_deaths2$plot.dat$Est)
# Linear model assumptions
z <- model.matrix(~ as.factor(year) * as.factor(month) + as.factor(dow), summer_data)
par(mfrow = c(1, 2))
plot(summer_deaths2$Yhat,
summer_data$logdeath - summer_deaths2$Yhat,
xlab = "Fitted", ylab = "Residuals"); abline(h = 0)
qqnorm(scale(summer_data$logdeath - summer_deaths2$Yhat, center = F), main = ""); abline(0, 1)
# Convergence
s <- lapply(mlist, summary, cenval = -20, pred.at = 10:32, mcmc = TRUE, verbose = FALSE)
summer_mcmc_tdlnm <- do.call(mcmc.list, lapply(1:restarts, function(i) {
d <- as.data.frame.table(s[[i]]$dlm_mcmc)
setDT(d)
d$Var1 <- paste0("Lag", rep(rep(0:lags, times = length(s[[i]]$pred.at)), times = iter/10))
d$Var2 <- paste0("Temp", rep(rep(s[[i]]$pred.at, each = length(0:lags)), times = iter/10))
mcmc(dcast(d[Var2 != "Temp20"], Var3 ~ Var1 + Var2, value.var = "Freq")[, Var3 := NULL][])
}))
gr_summer_tdlnm <- gelman.diag(summer_mcmc_tdlnm, transform = F, autoburnin = F, multivariate = F)GAM
cb <- crossbasis(summer_temp_dat, lag = c(0, lags),
argvar = list(fun = "ps", df = 10),
arglag = list(fun = "ps", intercept = T, df = 10))
pen <- cbPen(cb, addSlag = rep(0:1, c(3, 7)))
summer_deaths_gam <- gam(logdeath ~ cb + as.factor(year) * as.factor(month) + as.factor(dow),
data = summer_data, paraPen = list(cb = pen))
cp_summer <- crosspred(cb, summer_deaths_gam, at = 10:32, bylag = 1, cen = 20)
cp_summer$matfit <- out_trans(cp_summer$matfit)
cp_summer$matlow <- out_trans(cp_summer$matlow)
cp_summer$mathigh <- out_trans(cp_summer$mathigh)
# Linear model assumptions
par(mfrow = c(1, 2))
plot(summer_deaths_gam$fitted.values,
summer_data$logdeath - summer_deaths_gam$fitted.values,
xlab = "Fitted", ylab = "Residuals"); abline(h = 0)
qqnorm(scale(summer_deaths_gam$residuals, center = F), main = ""); abline(0, 1)
Lag days of susceptibility
We report periods of identified susceptiblity to temperature. For monotone-TDLNM, we use the posterior inference for probability of effect at each lag. For TDLNM and psDLNM we use periods where the credible intervals do not contain zero at any temperature level. This implies a relationship between temperature and mortality when contrasted with a temperature of 20 degrees C.
# Monotone TDLNM
plot(0:lags, s_summer_deaths$splitProb, main = "Monotone-TDLNM: Posterior inference", type = 'l', xlab = "Lag", ylab = "Prob of effect", ylim = c(0, 1)); abline(h = 0.95, col = "red", lty = 2)
(0:20)[s_summer_deaths$splitProb >= 0.95]## [1] 0 1 2 3# TDLNM
(0:20)[which(colSums(s_summer_deaths2$cilower > 0 | s_summer_deaths2$ciupper < 0) > 0)]## [1] 0 1 2 3 4 5 6 7 8 9 10 11 19# pSDLNM
(0:20)[which(colSums(cp_summer$matlow > 0 | cp_summer$mathigh < 0) > 0)]## [1] 0 1 2 3 4 5 6 8 9 10 11 13As before, we show the time periods and temperatures where there is a nonzero relationship between temperature and mortality. Blue indicates temperatures of increased mortality while red indicates temperatures with decreased mortality.
plot(s_summer_deaths, plot.type = "effect", main = "Monotone-TDLNM", start.time = 0, ylab = "Temperature (degrees C)", xlab = "Lag (days)")
plot(s_summer_deaths2, plot.type = "effect", main = "TDLNM", start.time = 0, ylab = "Temperature (degrees C)", xlab = "Lag (days)")
filled.contour(0:20, cp_summer$predvar, t((cp_summer$matlow>0) - (cp_summer$mathigh<0)),
xlab = "Lag (days)", ylab = "Temperature (degrees C)",
main = "GAM", levels = c(-1, -.1, .1, 1),
col = c("red", "white", "blue"))
Temperature slices of estimated exposure-lag-response function
plot_dat <- rbind.data.frame(
data.frame(lag = 0:lags,
y = s_summer_deaths$matfit["25",],
ymin = s_summer_deaths$cilower["25",],
ymax = s_summer_deaths$ciupper["25",],
grp = "25 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths$matfit["28",],
ymin = s_summer_deaths$cilower["28",],
ymax = s_summer_deaths$ciupper["28",],
grp = "28 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths$matfit["32",],
ymin = s_summer_deaths$cilower["32",],
ymax = s_summer_deaths$ciupper["32",],
grp = "32 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths2$matfit["25",],
ymin = s_summer_deaths2$cilower["25",],
ymax = s_summer_deaths2$ciupper["25",],
grp = "25 deg C",
model = "TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths2$matfit["28",],
ymin = s_summer_deaths2$cilower["28",],
ymax = s_summer_deaths2$ciupper["28",],
grp = "28 deg C",
model = "TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths2$matfit["32",],
ymin = s_summer_deaths2$cilower["32",],
ymax = s_summer_deaths2$ciupper["32",],
grp = "32 deg C",
model = "TDLNM"),
data.frame(lag = 0:lags,
y = cp_summer$matfit["25",],
ymin = cp_summer$matlow["25",],
ymax = cp_summer$mathigh["25",],
grp = "25 deg C",
model = "GAM"),
data.frame(lag = 0:lags,
y = cp_summer$matfit["28",],
ymin = cp_summer$matlow["28",],
ymax = cp_summer$mathigh["28",],
grp = "28 deg C",
model = "GAM"),
data.frame(lag = 0:lags,
y = cp_summer$matfit["32",],
ymin = cp_summer$matlow["32",],
ymax = cp_summer$mathigh["32",],
grp = "32 deg C",
model = "GAM")
)
plot_dat$model <- factor(plot_dat$model, levels = c("GAM", "TDLNM","Monotone-TDLNM"))
ggplot(plot_dat, aes(x = lag, y = y, ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 0, color = "red") +
geom_ribbon(fill = "grey") +
geom_line(size = 1) +
facet_grid(grp~model) +
theme_bw(base_size = 16) +
scale_x_continuous(expand = c(0, 0)) +
labs(x = "Lag (days)", y = "% Difference in Mortality")
Lag slices of estimated exposure-lag-response function
plot_dat <- rbind.data.frame(
data.frame(temp = s_summer_deaths$pred.at,
y = s_summer_deaths$matfit[,2],
ymin = s_summer_deaths$cilower[,2],
ymax = s_summer_deaths$ciupper[,2],
grp = "1 day prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_summer_deaths$pred.at,
y = s_summer_deaths$matfit[,3],
ymin = s_summer_deaths$cilower[,3],
ymax = s_summer_deaths$ciupper[,3],
grp = "2 days prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_summer_deaths$pred.at,
y = s_summer_deaths$matfit[,6],
ymin = s_summer_deaths$cilower[,6],
ymax = s_summer_deaths$ciupper[,6],
grp = "5 days prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_summer_deaths2$pred.at,
y = s_summer_deaths2$matfit[,2],
ymin = s_summer_deaths2$cilower[,2],
ymax = s_summer_deaths2$ciupper[,2],
grp = "1 day prior",
model = "TDLNM"),
data.frame(temp = s_summer_deaths2$pred.at,
y = s_summer_deaths2$matfit[,3],
ymin = s_summer_deaths2$cilower[,3],
ymax = s_summer_deaths2$ciupper[,3],
grp = "2 days prior",
model = "TDLNM"),
data.frame(temp = s_summer_deaths2$pred.at,
y = s_summer_deaths2$matfit[,6],
ymin = s_summer_deaths2$cilower[,6],
ymax = s_summer_deaths2$ciupper[,6],
grp = "5 days prior",
model = "TDLNM"),
data.frame(temp = cp_summer$predvar,
y = cp_summer$matfit[,2],
ymin = cp_summer$matlow[,2],
ymax = cp_summer$mathigh[,2],
grp = "1 day prior",
model = "GAM"),
data.frame(temp = cp_summer$predvar,
y = cp_summer$matfit[,3],
ymin = cp_summer$matlow[,3],
ymax = cp_summer$mathigh[,3],
grp = "2 days prior",
model = "GAM"),
data.frame(temp = cp_summer$predvar,
y = cp_summer$matfit[,6],
ymin = cp_summer$matlow[,6],
ymax = cp_summer$mathigh[,6],
grp = "5 days prior",
model = "GAM")
)
plot_dat$model <- factor(plot_dat$model, levels = c("GAM", "TDLNM","Monotone-TDLNM"))
ggplot(plot_dat, aes(x = temp, y = y, ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 0, color = "red") +
geom_ribbon(fill = "grey") +
geom_line(size = 1) +
facet_grid(grp~model) +
theme_bw(base_size = 16) +
scale_x_continuous(expand = c(0, 0)) +
labs(x = "Temperature", y = "% Difference in Mortality")
Heat maps
summer_plot_dat_gam <- s_summer_deaths$plot.dat
summer_plot_dat_gam$Est <- c(cp_summer$matfit)
summer_plot_dat <-
rbind.data.frame(cbind(s_summer_deaths$plot.dat, model = "Monotone-TDLNM"),
cbind(s_summer_deaths2$plot.dat, model = "TDLNM"),
cbind(summer_plot_dat_gam, model = "GAM"))
summer_plot_dat$model <- factor(summer_plot_dat$model, levels = c("GAM", "TDLNM","Monotone-TDLNM"))
ggplot(summer_plot_dat,
aes(xmin = Tmin - 1, xmax = Tmax - 1,
ymin = Xmin, ymax = Xmax, fill = Est)) +
geom_rect() +
scale_fill_viridis_c() +
scale_y_continuous(expand = c(0, 0)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw() +
facet_wrap(~model, nrow = 1) +
coord_equal(ratio = .65) +
labs(x = "Lag (days)", y = "Temperature", fill = "% Diff\nMortality", title = "")
Convergence results
Gelman-Rubin Statistic
# Effect
summer_rows <-
do.call(c, lapply(paste0("Lag", which(s_summer_deaths$splitProb > 0.95) - 1, "_"),
function(x) { grep(x, rownames(gr_summer_mtdlnm$psrf))}))
cv_summer <- data.frame("Monotone" = gr_summer_mtdlnm$psrf[summer_rows, 1],
"TDLNM" = gr_summer_tdlnm$psrf[summer_rows, 1])
cv_summer2 <- data.frame("Monotone" = gr_summer_mtdlnm$psrf[-summer_rows, 1],
"TDLNM" = gr_summer_tdlnm$psrf[-summer_rows, 1])
par(mfrow = c(1, 2))
boxplot(cv_summer, ylim = c(1, 2), main = "Summer: Effect", ylab = "Rhat"); abline(h = 1.2, lty = 2)
boxplot(cv_summer2, ylim = c(1, 2), main = "No Effect", ylab = "Rhat"); abline(h = 1.2, lty = 2)
Trace plots (Monotone-TDLNM)
summer_trace <- rbind(
rbindlist(lapply(1:restarts, function(i) {
rbindlist(lapply(c(15, 40, 50, 375), function(j) {
data.table(name = colnames(summer_mcmc_mtdlnm[[i]])[j],
model = "Monotone-TDLNM",
chain = i,
draws = 1:nrow(summer_mcmc_mtdlnm[[i]]),
sample = summer_mcmc_mtdlnm[[i]][, j])
}))
})),
rbindlist(lapply(1:restarts, function(i) {
rbindlist(lapply(c(15, 40, 50, 375), function(j) {
data.table(name = colnames(summer_mcmc_tdlnm[[i]])[j],
model = "TDLNM",
chain = i,
draws = 1:nrow(summer_mcmc_tdlnm[[i]]),
sample = summer_mcmc_tdlnm[[i]][, j])
}))
}))
)
ggplot(summer_trace) +
geom_line(aes(x = draws, y = sample, color = factor(chain)), alpha = 0.25) +
geom_density(aes(y = sample), orientation = "y",
alpha = 0.2, bw = 0.001) +
facet_grid(model~name) +
theme_minimal() +
theme(legend.position = "none") +
labs(x = "Draw", y = "Estimate")## Don't know how to automatically pick scale for object of type mcmc. Defaulting to continuous.
Winter model
We repeat the temperature-mortality model for winter days spanning October through April. We reverse the temperature axis to maintain monotone increasing for the monotone-TDLNM approach. All temperatures are represented as degrees below zero, where negative temperatures indicate degrees above zero. To incorporate less prior knowledge of the relationship between mortality and lagged cold temperatures, we specify a prior in Monotone-TDLNM indicating 95% confidence for the probability of an effect at any lag between 0.05 and 0.95.
Monotone-TDLNM
g0 <- rep(0, lags + 1) # prior prob .025-.975
s0 <- diag(lags + 1) * 2.701^2
ts0 <- rep(1/lags, lags)
mlist <-
foreach(m = 1:restarts, .errorhandling = "remove",
.packages = "dlmtree",
.verbose = FALSE) %dopar%
{
set.seed(m * 100)
tdlnm(logdeath ~ as.factor(year) * as.factor(month) + as.factor(dow),
data = winter_data,
exposure.data = winter_temp_dat,
exposure.splits = w_temp_splits,
n.trees = trees, n.burn = burn, n.iter = iter, n.thin = 10,
monotone.sigma = s0, altmcmc = 1,
monotone = T, verbose = F)
}
winter_deaths <- combine.models(mlist)
s_winter_deaths <- summary(winter_deaths, cenval = -20, pred.at = -20:20,
mcmc = TRUE, verbose = F)
s_winter_deaths$matfit <- out_trans(s_winter_deaths$matfit)
s_winter_deaths$cilower <- out_trans(s_winter_deaths$cilower)
s_winter_deaths$ciupper <- out_trans(s_winter_deaths$ciupper)
s_winter_deaths$plot.dat$Est <- out_trans(s_winter_deaths$plot.dat$Est)
par(mfrow = c(1, 2))
boxplot(1/(1+exp(-winter_deaths$zirtGamma)),
ylim = c(0, 1), main = "Prob. nonzero ERC", names = 0:lags);
points(1/(1+exp(-g0-1.96*sqrt(diag(s0)))), col = "red", pch = "-", cex = 3);
points(1/(1+exp(-g0+1.96*sqrt(diag(s0)))), col = "red", pch = "-", cex = 3);
boxplot(winter_deaths$timeProbs, ylim = c(0, 1),
main = "Prob. time split",
names = sapply(1:lags, function(i) paste0(i-1, "/", i)));
points(ts0, col = "blue", pch = 16)
# Linear model assumptions
par(mfrow = c(1, 2))
plot(winter_deaths$Yhat,
winter_data$logdeath - winter_deaths$Yhat,
xlab = "Fitted", ylab = "Residuals"); abline(h = 0)
qqnorm(scale(winter_data$logdeath - winter_deaths$Yhat, center = F), main = ""); abline(0, 1)
# Convergence
s <- lapply(mlist, summary, cenval = -20, pred.at = 10:32, mcmc = TRUE, verbose = FALSE)
winter_mcmc_mtdlnm <- do.call(mcmc.list, lapply(1:restarts, function(i) {
d <- as.data.frame.table(s[[i]]$dlm_mcmc)
setDT(d)
d$Var1 <- paste0("Lag", rep(rep(0:lags, times = length(s[[i]]$pred.at)), times = iter/10))
d$Var2 <- paste0("Temp", rep(rep(s[[i]]$pred.at, each = length(0:lags)), times = iter/10))
mcmc(dcast(d[Var2 != "Temp20"], Var3 ~ Var1 + Var2, value.var = "Freq")[, Var3 := NULL][])
}))
gr_winter_mtdlnm <- gelman.diag(winter_mcmc_mtdlnm, transform = F, autoburnin = F, multivariate = F)TDLNM
mlist <-
foreach(m = 1:restarts, .errorhandling = "remove",
.packages = "dlmtree",
.verbose = FALSE) %dopar%
{
set.seed(m * 100)
tdlnm(logdeath ~ as.factor(year) * as.factor(month) + as.factor(dow),
data = winter_data,
exposure.data = winter_temp_dat,
exposure.splits = w_temp_splits,
n.burn = burn, n.iter = iter, n.thin = 10,
verbose = F)
}
winter_deaths2 <- combine.models(mlist)
s_winter_deaths2 <- summary(winter_deaths2, cenval = -20, pred.at = -20:20,
mcmc = TRUE, verbose = F)
s_winter_deaths2$matfit <- out_trans(s_winter_deaths2$matfit)
s_winter_deaths2$cilower <- out_trans(s_winter_deaths2$cilower)
s_winter_deaths2$ciupper <- out_trans(s_winter_deaths2$ciupper)
s_winter_deaths2$plot.dat$Est <- out_trans(s_winter_deaths2$plot.dat$Est)
# Linear model assumptions
par(mfrow = c(1, 2))
plot(winter_deaths2$Yhat,
winter_data$logdeath - winter_deaths2$Yhat,
xlab = "Fitted", ylab = "Residuals"); abline(h = 0)
qqnorm(scale(winter_data$logdeath - winter_deaths2$Yhat, center = F), main = ""); abline(0, 1)
# Convergence
s <- lapply(mlist, summary, cenval = -20, pred.at = 10:32, mcmc = TRUE, verbose = FALSE)
winter_mcmc_tdlnm <- do.call(mcmc.list, lapply(1:restarts, function(i) {
d <- as.data.frame.table(s[[i]]$dlm_mcmc)
setDT(d)
d$Var1 <- paste0("Lag", rep(rep(0:lags, times = length(s[[i]]$pred.at)), times = iter/10))
d$Var2 <- paste0("Temp", rep(rep(s[[i]]$pred.at, each = length(0:lags)), times = iter/10))
mcmc(dcast(d[Var2 != "Temp20"], Var3 ~ Var1 + Var2, value.var = "Freq")[, Var3 := NULL][])
}))
gr_winter_tdlnm <- gelman.diag(winter_mcmc_tdlnm, transform = F, autoburnin = F, multivariate = F)GAM
cb <- crossbasis(winter_temp_dat, lag = c(0, lags),
argvar = list(fun = "ps", df = 10),
arglag = list(fun = "ps", intercept = T, df = 10))
pen <- cbPen(cb)
winter_deaths_gam <- gam(logdeath ~ cb + as.factor(year) * as.factor(month) + as.factor(dow),
data = winter_data, paraPen = list(cb = pen))
cp_winter <- crosspred(cb, winter_deaths_gam, at = -20:20, bylag = 1, cen = -20)
cp_winter$matfit <- out_trans(cp_winter$matfit)
cp_winter$matlow <- out_trans(cp_winter$matlow)
cp_winter$mathigh <- out_trans(cp_winter$mathigh)
# Linear model assumptions
par(mfrow = c(1, 2))
plot(winter_deaths_gam$fitted.values,
winter_data$logdeath - winter_deaths_gam$fitted.values,
xlab = "Fitted", ylab = "Residuals"); abline(h = 0)
qqnorm(scale(winter_deaths_gam$residuals, center = F), main = ""); abline(0, 1)
Lag days of susceptibility
# Monotone-TDLNM
plot(0:lags, s_winter_deaths$splitProb, main = "Monotone-TDLNM: Posterior inference", type = 'l', xlab = "Lag", ylab = "Prob of effect", ylim = c(0, 1)); abline(h = 0.95, col = "red", lty = 2)
(0:20)[s_winter_deaths$splitProb >= 0.95]## [1] 2 3 4 5 6 7 8 9# TDLNM
(0:20)[which(colSums(s_winter_deaths2$cilower > 0 | s_winter_deaths2$ciupper < 0) > 0)]## [1] 0 1 2 3 4 5 6 7 8 9 10# psDLNM
(0:20)[which(colSums(cp_winter$matlow > 0 | cp_winter$mathigh < 0) > 0)]## [1] 0 1 2 3 4 5 7 8 9plot(s_winter_deaths, plot.type = "effect", main = "Monotone-TDLNM", start.time = 0, ylab = "Temperature (degrees below C)", xlab = "Lag (days)")
plot(s_winter_deaths2, plot.type = "effect", main = "TDLNM", start.time = 0, ylab = "Temperature (degrees below C)", xlab = "Lag (days)")
filled.contour(0:20, cp_winter$predvar, t((cp_winter$matlow>0) - (cp_winter$mathigh<0)),
xlab = "Lag (days)", ylab = "Temperature (degrees below C)",
main = "GAM", levels = c(-1, -.1, .1, 1),
col = c("red", "white", "blue"))
Temperature slices of estimated exposure-lag-response function
plot_dat <- rbind.data.frame(
data.frame(lag = 0:lags,
y = s_winter_deaths$matfit["-10",],
ymin = s_winter_deaths$cilower["-10",],
ymax = s_winter_deaths$ciupper["-10",],
grp = "10 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_winter_deaths$matfit["0",],
ymin = s_winter_deaths$cilower["0",],
ymax = s_winter_deaths$ciupper["0",],
grp = "0 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_winter_deaths$matfit["15",],
ymin = s_winter_deaths$cilower["15",],
ymax = s_winter_deaths$ciupper["15",],
grp = "-15 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_winter_deaths2$matfit["-10",],
ymin = s_winter_deaths2$cilower["-10",],
ymax = s_winter_deaths2$ciupper["-10",],
grp = "10 deg C",
model = "TDLNM"),
data.frame(lag = 0:lags,
y = s_winter_deaths2$matfit["0",],
ymin = s_winter_deaths2$cilower["0",],
ymax = s_winter_deaths2$ciupper["0",],
grp = "0 deg C",
model = "TDLNM"),
data.frame(lag = 0:lags,
y = s_winter_deaths2$matfit["15",],
ymin = s_winter_deaths2$cilower["15",],
ymax = s_winter_deaths2$ciupper["15",],
grp = "-15 deg C",
model = "TDLNM"),
data.frame(lag = 0:lags,
y = cp_winter$matfit["-10",],
ymin = cp_winter$matlow["-10",],
ymax = cp_winter$mathigh["-10",],
grp = "10 deg C",
model = "GAM"),
data.frame(lag = 0:lags,
y = cp_winter$matfit["0",],
ymin = cp_winter$matlow["0",],
ymax = cp_winter$mathigh["0",],
grp = "0 deg C",
model = "GAM"),
data.frame(lag = 0:lags,
y = cp_winter$matfit["15",],
ymin = cp_winter$matlow["15",],
ymax = cp_winter$mathigh["15",],
grp = "-15 deg C",
model = "GAM")
)
plot_dat$model <- factor(plot_dat$model, levels = c("GAM", "TDLNM","Monotone-TDLNM"))
ggplot(plot_dat, aes(x = lag, y = y, ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 0, color = "red") +
geom_ribbon(fill = "grey") +
geom_line(size = 1) +
facet_grid(grp~model) +
theme_bw(base_size = 16) +
scale_x_continuous(expand = c(0, 0)) +
labs(x = "Lag (days)", y = "% Difference in Mortality")
Lag slices of estimated exposure-lag-response function
plot_dat <- rbind.data.frame(
data.frame(temp = s_winter_deaths$pred.at,
y = s_winter_deaths$matfit[,2],
ymin = s_winter_deaths$cilower[,2],
ymax = s_winter_deaths$ciupper[,2],
grp = "1 day prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_winter_deaths$pred.at,
y = s_winter_deaths$matfit[,3],
ymin = s_winter_deaths$cilower[,3],
ymax = s_winter_deaths$ciupper[,3],
grp = "2 days prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_winter_deaths$pred.at,
y = s_winter_deaths$matfit[,6],
ymin = s_winter_deaths$cilower[,6],
ymax = s_winter_deaths$ciupper[,6],
grp = "5 days prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_winter_deaths2$pred.at,
y = s_winter_deaths2$matfit[,2],
ymin = s_winter_deaths2$cilower[,2],
ymax = s_winter_deaths2$ciupper[,2],
grp = "1 day prior",
model = "TDLNM"),
data.frame(temp = s_winter_deaths2$pred.at,
y = s_winter_deaths2$matfit[,3],
ymin = s_winter_deaths2$cilower[,3],
ymax = s_winter_deaths2$ciupper[,3],
grp = "2 days prior",
model = "TDLNM"),
data.frame(temp = s_winter_deaths2$pred.at,
y = s_winter_deaths2$matfit[,6],
ymin = s_winter_deaths2$cilower[,6],
ymax = s_winter_deaths2$ciupper[,6],
grp = "5 days prior",
model = "TDLNM"),
data.frame(temp = cp_winter$predvar,
y = cp_winter$matfit[,2],
ymin = cp_winter$matlow[,2],
ymax = cp_winter$mathigh[,2],
grp = "1 day prior",
model = "GAM"),
data.frame(temp = cp_winter$predvar,
y = cp_winter$matfit[,3],
ymin = cp_winter$matlow[,3],
ymax = cp_winter$mathigh[,3],
grp = "2 days prior",
model = "GAM"),
data.frame(temp = cp_winter$predvar,
y = cp_winter$matfit[,6],
ymin = cp_winter$matlow[,6],
ymax = cp_winter$mathigh[,6],
grp = "5 days prior",
model = "GAM")
)
plot_dat$model <- factor(plot_dat$model, levels = c("GAM", "TDLNM","Monotone-TDLNM"))
ggplot(plot_dat, aes(x = temp, y = y, ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 0, color = "red") +
geom_ribbon(fill = "grey") +
geom_line(size = 1) +
facet_grid(grp~model) +
theme_bw(base_size = 16) +
scale_x_continuous(expand = c(0, 0), limits = c(-19, 20)) +
labs(x = "Degrees below zero", y = "% Difference in Mortality")## Warning: Removed 1 row(s) containing missing values (geom_path).
Heat maps
winter_plot_dat_gam <- s_winter_deaths$plot.dat
winter_plot_dat_gam$Est <- c(cp_winter$matfit)
winter_plot_dat <-
rbind.data.frame(cbind(s_winter_deaths$plot.dat, model = "Monotone-TDLNM"),
cbind(s_winter_deaths2$plot.dat, model = "TDLNM"),
cbind(winter_plot_dat_gam, model = "GAM"))
winter_plot_dat$model <- factor(winter_plot_dat$model, levels = c("GAM", "TDLNM","Monotone-TDLNM"))
ggplot(winter_plot_dat,
aes(xmin = Tmin - 1, xmax = Tmax - 1,
ymin = Xmin, ymax = Xmax, fill = Est)) +
geom_rect() +
scale_fill_viridis_c() +
scale_y_continuous(expand = c(0, 0)) + scale_x_continuous(expand = c(0, 0)) +
theme_bw() +
facet_wrap(~model, nrow = 1) +
coord_equal(ratio = .4) +
labs(x = "Lag (days)", y = "Degrees below zero", fill = "% Diff\nMortality", title = "")
Convergence results
Gelman-Rubin Statistic
# Effect
winter_rows <-
do.call(c, lapply(paste0("Lag", which(s_winter_deaths$splitProb > 0.9) - 1, "_"),
function(x) { grep(x, rownames(gr_winter_mtdlnm$psrf))}))
cv_winter <- data.frame("Monotone" = gr_winter_mtdlnm$psrf[winter_rows, 1],
"TDLNM" = gr_winter_tdlnm$psrf[winter_rows, 1])
cv_winter2 <- data.frame("Monotone" = gr_winter_mtdlnm$psrf[-winter_rows, 1],
"TDLNM" = gr_winter_tdlnm$psrf[-winter_rows, 1])
par(mfrow = c(1, 2))
boxplot(cv_winter, ylim = c(1, 2), main = "Winter: Effect", ylab = "Rhat"); abline(h = 1.2, lty = 2)
boxplot(cv_winter2, ylim = c(1, 2), main = "No Effect", ylab = "Rhat"); abline(h = 1.2, lty = 2)
Trace plots (Monotone-TDLNM)
winter_trace <- rbind(
rbindlist(lapply(1:restarts, function(i) {
rbindlist(lapply(c(20, 230, 275, 175), function(j) {
data.table(name = colnames(winter_mcmc_mtdlnm[[i]])[j],
model = "Monotone-TDLNM",
chain = i,
draws = 1:nrow(winter_mcmc_mtdlnm[[i]]),
sample = winter_mcmc_mtdlnm[[i]][, j])
}))
})),
rbindlist(lapply(1:restarts, function(i) {
rbindlist(lapply(c(20, 230, 275, 175), function(j) {
data.table(name = colnames(winter_mcmc_tdlnm[[i]])[j],
model = "TDLNM",
chain = i,
draws = 1:nrow(winter_mcmc_tdlnm[[i]]),
sample = winter_mcmc_tdlnm[[i]][, j])
}))
}))
)
ggplot(winter_trace) +
geom_line(aes(x = draws, y = sample, color = factor(chain)), alpha = 0.25) +
geom_density(aes(y = sample), orientation = "y",
alpha = 0.2, bw = 0.001) +
facet_grid(model~name) +
theme_minimal() +
theme(legend.position = "none") +
labs(x = "Draw", y = "Estimate")## Don't know how to automatically pick scale for object of type mcmc. Defaulting to continuous.
Probability of effect
plot_dat_pe <- rbind.data.frame(
data.frame(lag = 0:lags, prob = s_summer_deaths$splitProb, group = "Summer"),
data.frame(lag = 0:lags, prob = s_winter_deaths$splitProb, group = "Winter")
)
ggplot(plot_dat_pe, aes(x = lag, y = prob, color = group, linetype = group)) +
geom_hline(yintercept = 0.95, color = "black", linetype = 1, size = 1) +
geom_line(size = 2) +
theme_bw(base_size = 16) +
scale_y_continuous(limits = c(0, 1)) +
scale_x_continuous(expand = c(0, 0)) +
labs(x = "Lag (days)", y = "P(susceptibility at lag l)", color = "", linetype = "") +
theme(legend.position = "bottom", legend.key.width = unit(2, "cm"))
Convergence results
cv_summer <- data.frame("Monotone" = gr_summer_mtdlnm$psrf[, 1],
"TDLNM" = gr_summer_tdlnm$psrf[, 1])
cv_winter <- data.frame("Monotone" = gr_winter_mtdlnm$psrf[, 1],
"TDLNM" = gr_winter_tdlnm$psrf[, 1])
par(mfrow = c(1, 2))
boxplot(cv_summer, ylim = c(1, 2), main = "Gelman-Rubin Statistic:\nSummer", ylab = "Rhat"); abline(h = 1.2, lty = 2)
boxplot(cv_winter2, ylim = c(1, 2), main = "\nWinter", ylab = "Rhat"); abline(h = 1.2, lty = 2)
Sensitivity analyses
In a sensitivity analysis for our summer mortality dataset, we rerun Monotone-TDLNM but with the vague prior used in the winter mortality analysis.
Summer: Monotone-TDLNM, vague prior
s0 <- diag(lags + 1) * 2.701^2
mlist <-
foreach(m = 1:restarts, .errorhandling = "remove",
.packages = "dlmtree",
.verbose = FALSE) %dopar%
{
set.seed(m * 100)
tdlnm(logdeath ~ as.factor(year) * as.factor(month) + as.factor(dow),
data = summer_data,
exposure.data = summer_temp_dat,
exposure.splits = temp_splits,
n.trees = trees, n.burn = burn, n.iter = iter, n.thin = 10,
monotone.sigma = s0, altmcmc = 1,
monotone = T, verbose = F)
}
stopCluster(cl)
summer_deaths0 <- combine.models(mlist)
s_summer_deaths0 <- summary(summer_deaths0, cenval = 20, pred.at = 10:32,
mcmc = TRUE, verbose = F)
s_summer_deaths0$matfit <- out_trans(s_summer_deaths0$matfit)
s_summer_deaths0$cilower <- out_trans(s_summer_deaths0$cilower)
s_summer_deaths0$ciupper <- out_trans(s_summer_deaths0$ciupper)
s_summer_deaths0$plot.dat$Est <- out_trans(s_summer_deaths0$plot.dat$Est)plot(0:lags, s_summer_deaths0$splitProb, main = "Monotone-TDLNM: Posterior inference", type = 'l', xlab = "Lag", ylab = "Prob of effect", ylim = c(0, 1)); abline(h = 0.95, col = "red", lty = 2)
plot_dat <- rbind.data.frame(
data.frame(lag = 0:lags,
y = s_summer_deaths$matfit["25",],
ymin = s_summer_deaths$cilower["25",],
ymax = s_summer_deaths$ciupper["25",],
grp = "25 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths$matfit["28",],
ymin = s_summer_deaths$cilower["28",],
ymax = s_summer_deaths$ciupper["28",],
grp = "28 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths$matfit["32",],
ymin = s_summer_deaths$cilower["32",],
ymax = s_summer_deaths$ciupper["32",],
grp = "32 deg C",
model = "Monotone-TDLNM"),
data.frame(lag = 0:lags,
y = s_summer_deaths0$matfit["25",],
ymin = s_summer_deaths0$cilower["25",],
ymax = s_summer_deaths0$ciupper["25",],
grp = "25 deg C",
model = "Uninformative prior"),
data.frame(lag = 0:lags,
y = s_summer_deaths0$matfit["28",],
ymin = s_summer_deaths0$cilower["28",],
ymax = s_summer_deaths0$ciupper["28",],
grp = "28 deg C",
model = "Uninformative prior"),
data.frame(lag = 0:lags,
y = s_summer_deaths0$matfit["32",],
ymin = s_summer_deaths0$cilower["32",],
ymax = s_summer_deaths0$ciupper["32",],
grp = "32 deg C",
model = "Uninformative prior")
)
ggplot(plot_dat, aes(x = lag, y = y, ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 0, color = "red") +
geom_ribbon(fill = "grey") +
geom_line(size = 1) +
facet_grid(model ~ grp) +
theme_bw(base_size = 16) +
scale_x_continuous(expand = c(0, 0)) +
labs(x = "Lag (days)", y = "% Difference in Mortality")
plot_dat <- rbind.data.frame(
data.frame(temp = s_summer_deaths$pred.at,
y = s_summer_deaths$matfit[,2],
ymin = s_summer_deaths$cilower[,2],
ymax = s_summer_deaths$ciupper[,2],
grp = "1 day prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_summer_deaths$pred.at,
y = s_summer_deaths$matfit[,3],
ymin = s_summer_deaths$cilower[,3],
ymax = s_summer_deaths$ciupper[,3],
grp = "2 days prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_summer_deaths$pred.at,
y = s_summer_deaths$matfit[,6],
ymin = s_summer_deaths$cilower[,6],
ymax = s_summer_deaths$ciupper[,6],
grp = "5 days prior",
model = "Monotone-TDLNM"),
data.frame(temp = s_summer_deaths0$pred.at,
y = s_summer_deaths0$matfit[,2],
ymin = s_summer_deaths0$cilower[,2],
ymax = s_summer_deaths0$ciupper[,2],
grp = "1 day prior",
model = "Uninformative prior"),
data.frame(temp = s_summer_deaths0$pred.at,
y = s_summer_deaths0$matfit[,3],
ymin = s_summer_deaths0$cilower[,3],
ymax = s_summer_deaths0$ciupper[,3],
grp = "2 days prior",
model = "Uninformative prior"),
data.frame(temp = s_summer_deaths0$pred.at,
y = s_summer_deaths0$matfit[,6],
ymin = s_summer_deaths0$cilower[,6],
ymax = s_summer_deaths0$ciupper[,6],
grp = "5 days prior",
model = "Uninformative prior")
)
ggplot(plot_dat, aes(x = temp, y = y, ymin = ymin, ymax = ymax)) +
geom_hline(yintercept = 0, color = "red") +
geom_ribbon(fill = "grey") +
geom_line(size = 1) +
facet_grid(model ~ grp) +
theme_bw(base_size = 16) +
scale_x_continuous(expand = c(0, 0)) +
labs(x = "Temperature", y = "% Difference in Mortality")